﻿using LxBasic;
using LxLESS;

namespace LxCore
{
    public partial class MatrixBase : ScalarBase, IZoneIteration
    {
        public double[] ADiag;//线性方程组对角线
        public double[] B;//线性方程组 b
        public double[] AInnO;
        public double[] AInnN;
        public double[] ABouH;//FormAB阶段用于添加到ADiag的计算中，然后立即清零，矩阵求解中用于不同zone交换数据


        public IParameter<double> RelaxFactor = ConstField.Double1;

        #region 初始化及公用
        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);
            CheckCreator.ArrayStrict(ref B, InnCells.Length);
            CheckCreator.ArrayStrict(ref ADiag, InnCells.Length);
            CheckCreator.ArrayStrict(ref AInnO, InnFaces.Length);
            if (IsSymmetric())
            {
                AInnN = AInnO;
            }
            else
            {
                if (ReferenceEquals(AInnO, AInnN)) { AInnN = null; }
                CheckCreator.ArrayStrict(ref AInnN, InnFaces.Length);
            }
            CheckCreator.ArrayStrict(ref ABouH, BouFaces.Length);
        }
        public virtual bool IsSymmetric()
        {
            return false;
        }
        public override void Reuse(ZoneHead oldOne)
        {
            base.Reuse(oldOne);
            if (oldOne is MatrixBase o)
            {
                this.ADiag = o.ADiag;
                this.B = o.B;
                this.AInnO = o.AInnO;
                this.AInnN = o.AInnN;
            }
        }
        public virtual void FormABOnFaces(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
        }
        public virtual void FormABOnCells(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
        }
        #endregion
        #region 计算调用
        protected unsafe double CalAdiag(ref FvInnCell cell)
        {
            double result = 0;
            fixed (int* start = cell.AdjFaces)
            {
                int* adjf = start, endp1 = start + cell.InnOFaceEnd;
                for (; adjf < endp1; ++adjf)
                {
                    result += AInnN[*adjf];//加对面的
                }
                endp1 = start + cell.InnNFaceEnd;
                for (; adjf < endp1; ++adjf)
                {
                    result += AInnO[*adjf];
                }
                endp1 = start + cell.AdjFaces.Length;
                for (; adjf < endp1; ++adjf)
                {
                    result += ABouH[*adjf];
                    ABouH[*adjf] = 0.0;//同时清空，用于后面存储区域间传递的数据
                }
            }
            return -result;//为负值
        }
        protected void FormALUb_0(IIndexHost indexNumber, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = indexNumber.FaceIdx.Split(threadInfo);
            int len = endp1 - i1;
            if (indexNumber is FvInner)
            {
                Array.Clear(AInnO, i1, len);
                Array.Clear(AInnN, i1, len);
                Array.Clear(InnFaceShare, i1, len);
            }
            else if (indexNumber is FvBound)
            {
                Array.Clear(ABouH, i1, len);
                Array.Clear(BouFaceShare, i1, len);
            }
        }
        protected unsafe double calALUProductValue(ref FvInnCell cell, double[] fromArray)
        {
            double aluProductValue = 0;
            fixed (int* start = cell.AdjFaces)
            {
                int* adjf = start, endp1 = start + cell.InnOFaceEnd;
                for (; adjf < endp1; ++adjf)
                {
                    aluProductValue += AInnO[*adjf] * fromArray[InnFaces[*adjf].NSideCell];
                }
                endp1 = start + cell.InnNFaceEnd;
                for (; adjf < endp1; ++adjf)
                {
                    aluProductValue += AInnN[*adjf] * fromArray[InnFaces[*adjf].OSideCell];
                }
                endp1 = start + cell.AdjFaces.Length;
                for (; adjf < endp1; ++adjf)
                {
                    aluProductValue += ABouH[*adjf];
                }
            }
            return aluProductValue;
        }
        protected void cal_x(double[] fromArray, double[] toArray, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                toArray[i1] = (B[i1] - calALUProductValue(ref InnCells[i1], fromArray)) / ADiag[i1];
            }
        }
        #endregion
        #region Relax
        public void ImplicitUnderRelax(ThreadInfo threadInfo)
        {
            ImplicitUnderRelax(RelaxFactor.GetPara(), threadInfo);
        }
        public void ImplicitUnderRelax(double relaxFactor, ThreadInfo threadInfo)
        {
            if (relaxFactor == 1) return;
            double relax1 = 1 - relaxFactor;
            double relaxFactorReverse = 1 / relaxFactor;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref double adiag = ref ADiag[i1];
                adiag *= relaxFactorReverse;
                B[i1] += relax1 * adiag * InnX[i1];
            }
        }
        //private void ImplicitUnderRelax_bPart(ThreadInfo threadInfo)//当A不变时，relax不需要重复求A，仅求b减少计算量
        //{
        //    var rf = RelaxFactor.GetCellVars(null);
        //    double relaxFactor = rf.GetCellVar(LxSim.Sim.Steps.SteadyStep.StepsCount);
        //    double relax1 = (1 - relaxFactor);
        //    (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
        //    for (; i1 < endp1; i1++)
        //    {
        //        B[i1] += relax1 * ADiag[i1] * InnX[i1];
        //    }
        //}

        #endregion
        #region IZoneIteration
        protected RMSResidualTool Residual = new RMSResidualTool();
        public void CalResidual(ThreadInfo threadInfo)
        {
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            RMSResidualTool res = new RMSResidualTool();
            for (; i1 < endp1; i1++)
            {
                double future_x = (B[i1] - calALUProductValue(ref InnCells[i1], InnX)) / ADiag[i1];
                res.Include(InnX[i1], future_x);
            }
            lock (MultiThreads.GlobalLock)
            {
                Residual.Include(res);
            }
        }

        public (bool isValid, double res) GetResidualReset()
        {
            double res = Residual.GetRMS(InnCells.Length);
            Residual = new RMSResidualTool();
            return (true, res);
        }

        public void JacobiIterate_Step1(ThreadInfo threadInfo)
        {
            cal_x(InnX, InnFaceShare, threadInfo);
        }
        public void JacobiIterate_Step2(ThreadInfo threadInfo)
        {
            cal_x(InnFaceShare, InnX, threadInfo);
        }
        public void ExplicitIterate_Step1(ThreadInfo threadInfo)
        {
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            Array.Copy(InnX, i1, InnFaceShare, i1, endp1 - i1);
        }
        public void ExplicitIterate_Step2(ThreadInfo threadInfo)
        {
            cal_x(InnFaceShare, InnX, threadInfo);
        }
        public void AMGIte_Restrict_Are(ThreadInfo threadInfo)
        {
            MGZone coarse = GetCoarserZone();// Zone.AMGSolver.MgZones[0];
            var coarseCells = coarse.InnCells;
            var coarsee = coarse.e;
            (int i1, int endp1) = coarse.InnerCellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref MGCell coarCell = ref coarseCells[i1];
                coarCell.ADiag = restrict_Ar(ref coarCell);
                coarsee[i1] = coarCell.r / coarCell.ADiag;
            }
            //A_LU
            coarse.HelpRestrict_ALU(AInnO, AInnN, threadInfo);
        }

        public void AMGIte_Restrict_re(ThreadInfo threadInfo)
        {
            MGZone coarse = GetCoarserZone();// Zone.AMGSolver.MgZones[0];
            var coarseCells = coarse.InnCells;
            var coarsee = coarse.e;
            (int i1, int endp1) = coarse.InnerCellIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref MGCell coarCell = ref coarseCells[i1];
                restrict_r(ref coarCell);
                coarsee[i1] = coarCell.r / coarCell.ADiag;
            }
        }
        public void AMGIte_ReceiveProlong(ThreadInfo threadInfo)
        {
            GetCoarserZone().HelpProlong(InnX, threadInfo);
        }
        private double restrict_Ar(ref MGCell mgCell)
        {
            double Atemp = 0;
            double rtemp = 0;
            foreach (int finecell in mgCell.FineCells)
            {
                double atemp = ADiag[finecell];
                Atemp += atemp;
                rtemp += B[finecell] - calALUProductValue(ref InnCells[finecell], InnX) - atemp * InnX[finecell];
            }
            foreach (int innerfineface in mgCell.InnerFineFaces)
            {
                Atemp += AInnO[innerfineface] + AInnN[innerfineface];
            }
            mgCell.r = rtemp;
            return Atemp;
        }
        private void restrict_r(ref MGCell coarseCell)
        {
            double rtemp = 0;
            foreach (int finecell in coarseCell.FineCells)
            {
                rtemp += B[finecell] - calALUProductValue(ref InnCells[finecell], InnX) - ADiag[finecell] * InnX[finecell];
            }
            coarseCell.r = rtemp;
        }

        public MGZone GetCoarserZone()
        {
            return GetLessZone().MGZones[0];
        }
        public LESSZone GetLessZone()
        {
            return Zone.LESSZone;
        }
        public void BiCGStabIte_Ini_r_rho1(ThreadInfo threadInfo)
        {
            var r = GetBiCGStabSolver().r;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            Float128 temp = new();
            for (; i1 < endp1; i1++)
            {
                double rt = r[i1] = B[i1] - calALUProductValue(ref InnCells[i1], InnX) - ADiag[i1] * InnX[i1];
                temp += (rt * rt);
            }
            ThreadsMath.MemoryAdd(temp);
        }
        public void BiCGStabIte_Ini_r0_pArray_yArray(ThreadInfo threadInfo)
        {
            GetBiCGStabSolver().BiCGStab_JacPrec_Ini_r0_pArray_yArray(Zone.Inner.CellIdx, ADiag, threadInfo);
        }
        public void BiCGStabIte_vArray_alpha1(ThreadInfo threadInfo)
        {
            BiCGStabZone bicg = GetBiCGStabSolver();
            double[] v = bicg.v, y = bicg.y, r0 = bicg.r0;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            Float128 temp = new();
            for (; i1 < endp1; i1++)
            {
                double vt = v[i1] = ADiag[i1] * y[i1] + calALUProductValue(ref InnCells[i1], y);
                temp += r0[i1] * vt;
            }
            ThreadsMath.MemoryAdd(temp);
        }
        public void BiCGStabIte_sArray_zArray(ThreadInfo threadInfo)
        {
            GetBiCGStabSolver().BiCGStab_JacPrec_sArray_zArray(Zone.Inner.CellIdx, ADiag, threadInfo);
        }
        public void BiCGStabIte_tArray_omega1(ThreadInfo threadInfo)
        {
            BiCGStabZone bicg = GetBiCGStabSolver();
            double[] s = bicg.r, z = bicg.z, t = bicg.t;
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            Float128 temp1 = new(), temp2 = new();
            for (; i1 < endp1; i1++)
            {
                double tt = t[i1] = ADiag[i1] * z[i1] + calALUProductValue(ref InnCells[i1], z);
                temp1 += tt * s[i1];
                temp2 += tt * tt;
            }
            ThreadsMath.MemoryAdd(temp1, temp2);
        }
        public void BiCGStabIte_xArray(ThreadInfo threadInfo)
        {
            GetBiCGStabSolver().BiCGStab_Cal_xArray(Zone.Inner.CellIdx, InnX, threadInfo);
        }
        public void BiCGStabIte_rArray_rho1(ThreadInfo threadInfo)
        {
            GetBiCGStabSolver().BiCGStab_rArray_rho1(Zone.Inner.CellIdx, threadInfo);
        }
        public void BiCGStabIte_pArray_yArray(ThreadInfo threadInfo)
        {
            GetBiCGStabSolver().BiCGStab_JacPrec_pArray_yArray(Zone.Inner.CellIdx, ADiag, threadInfo);
        }
        public BiCGStabZone GetBiCGStabSolver()
        {
            return Zone.LESSZone.BiCGZone;
        }
        #endregion
    }
}
//public class MatrixBoundBase : MatrixBoundBase
//{
//    //public override void Initialize(object inMatrix, FvBound fvThread, ThreadInfo threadInfo)
//    //{
//    //    var m = inMatrix as MatrixBase;
//    //    m.InitializeX(m.IniValue, fvThread, threadInfo);
//    //}

//}  //public virtual void Initialize(object inMatrix, FvBound fvBound, ThreadInfo threadInfo) { }
//public virtual void FormBoundValue(object inMatrix, FvBound fvBound, ThreadInfo threadInfo) { }
//public virtual void FormAbOnFaces(object inMatrix, FvBound fvBound, ThreadInfo threadInfo) { }
//public virtual void UpdateBoundValue(object inMatrix, FvBound fvBound, ThreadInfo threadInfo) { }